

#include "slep.h"
#include "hr_time.h"
#include "epph.h"
#include "float.h"

FUNVAL mtLeastRm1(double* A, double* y, double z, OPTS opt);
FUNVAL mtLeastRm2(double* A, double* y, double z, OPTS opt);
FUNVAL mtLeastRm3(double* A, double* y, double z, OPTS opt);


FUNVAL mtLeastR(double* A, double* y, double z, OPTS opt){
	if(opt.mFlag== 0 && opt.lFlag == 0 )
		return mtLeastRm1(A,y,z,opt);
	else if(opt.mFlag == 1 && opt.lFlag == 0 && opt.q == 2)
		return mtLeastRm2(A,y,z,opt);
	else if(opt.mFlag == 1 && opt.lFlag == 1 && opt.q == 2)
		return mtLeastRm3(A,y,z,opt);
	else
		errorOut("Please set the values of opt.mFlag and opt.lFlag to correct values before executing mtLeastR.");
}


FUNVAL mtLeastRm1(double* A, double* y,double z,OPTS opt){

	FUNVAL ret;
	int m,n,k,i, iterStep;
	int *ind;
	double q, beta, alpha, alphap, q_bar, lambda, lambda_max, L, l_sum, r_sum;
	double fVal, fValp;
	double *tt, *ATy,*x, *Ax, *xxp, *xp, *s, *As, *ATAs, *g, *Axp, *v, *Av;
	double *ValueL, *funVal, *Axy;

	m = opt.m;
	n = opt.n;
	if(!hasGot(&opt,"ind"))
		errorOut("In mtLeastR, 'ind' should be specified.");
	else{
		ind = opt.ind;
		k = opt.ind_size - 1;
		if(ind[k] != (m))
			errorOut("Check opt.ind and its size.");
	}

	tt = (double *) malloc( sizeof(double) * n);
	ATy = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc( sizeof(double) * n*k);
	Ax = (double *) malloc( sizeof(double) * m);
	xp = (double *) malloc( sizeof(double) * n*k);
	Axp = (double *) malloc( sizeof(double) * m);
	xxp = (double *) malloc( sizeof(double) * n*k);
	s = (double *) malloc( sizeof(double) * n*k);
	As = (double *) malloc( sizeof(double) * m);
	ATAs = (double *) malloc( sizeof(double) * n*k);
	g = (double *) malloc( sizeof(double) * n*k);
	v = (double *) malloc( sizeof(double) * n*k);
	Av = (double *) malloc( sizeof(double) * m);
	ValueL = (double *) malloc( sizeof(double) * opt.maxIter);
	Axy = (double *) malloc( sizeof(double) * m);
	funVal = (double *) malloc( sizeof(double) * opt.maxIter);

	if(!hasGot(&opt,"q"))
		opt.q = q = 2;
	else if( (q = opt.q) < 1)
		errorOut("'q' should be larger than 1.");
	
	initNormalization(&opt,A);

	/******* Starting point initialization *******/
/*
ATy=zeros(n, k);
% compute AT y
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    tt =A(ind_i,:)'*y(ind_i,1); / normalized
    ATy(:,i)= tt;
end
*/
	for(i=0; i < k; i++){
		normalizedmv_lowlevel('T',ind[i+1]-ind[i] ,n,&A[ind[i]],m,&y[ind[i]],tt,ind[i],0, &opt);
		dcopy(n,tt,1,&ATy[i*n],1);
	}

	// process the regularization parameter
	if(opt.rFlag == 0)
		lambda = z;
	else{
		if( z<0.0 || z>1.0)
			errorOut("Z is not in range (0,1)");
		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATy[i],n,q_bar));
		lambda = z * lambda_max;
	}

	if( opt.init == 2)
		memset(x,0,sizeof(double)*n*k);
	else{
		if( hasGot(&opt,"x0"))
			x = opt.x0;
		else
			x = ATy;
	}

	memset(Ax,0,sizeof(double)*m);
	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);

	if(opt.init == 0) {
		// TO-DO
	}

	/*........ Method specific code starts here .........*/

	/************************ Armijo Line Search *****/
	if (opt.mFlag==0 && opt.lFlag==0){
		L=1.0;

		//xp=x; Axp=Ax; xxp=zeros(n,k);
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		memset(xxp,0,n*k*sizeof(double));

		alphap=0.0; alpha=1.0;

		for(iterStep = 0; iterStep < opt.maxIter ; iterStep ++ ){
			beta=(alphap-1)/alpha;    
			//s=x + beta* xxp;
			dcopy(n*k,x,1,s,1);
			daxpy(n*k,beta,xxp,1,s,1);

			//As=Ax + beta* (Ax-Axp);
			dcopy(m,Ax,1,As,1);
			dscal(m,1.0 + beta, As,1);
			daxpy(m,-beta,Axp,1,As,1);

			for( i = 0 ; i<k ; i++){
				normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&As[ind[i]],tt,ind[i],0, &opt);
				dcopy(n,tt,1,&ATAs[i*n],1);
			}
			//g=ATAs-ATy;
			dcopy(n*k,ATAs,1,g,1);
			daxpy(n*k,-1.0,ATy,1,g,1);

			//xp=x;    Axp=Ax;
			dcopy(n*k,x,1,xp,1);
			dcopy(m,Ax,1,Axp,1);

			for(;;){
				//v=s-g/L;
				dcopy(n*k,s,1,v,1);
				daxpy(n*k,-1.0/L,g,1,v,1);

				//x=eppMatrix(v, n, k, lambda/ L, q);
				//eppMatrix(x,v,k,n,lambda/L,q);
				eppMatrix(x,v,n,k,lambda/L,q);

				// v=x-s;
				dcopy(n*k,x,1,v,1);
				daxpy(n*k,-1.0,s,1,v,1);

				memset(Ax,0,sizeof(double)*m);
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0, &opt);
				
				//Av=Ax -As;
				dcopy(m,Ax,1,Av,1);
				daxpy(m,-1.0,As,1,Av,1);

				//r_sum=norm(v,'fro')^2; l_sum=Av'*Av;
				r_sum = ddot(n*k,v,1,v,1);
				l_sum = ddot(m,Av,1,Av,1);

				if(l_sum <= r_sum * L)
					break;
				else
					L = maxof2(2.0*L, l_sum/r_sum);
			}

			alphap=alpha; 
			alpha= (1.0l+ sqrt(4.0l*alpha*alpha +1.0l))/2.0l;

			ValueL [iterStep] = L;
			// xxp = x - xp;
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1.0,xp,1,xxp,1);

			//Axy=Ax-y;
			dcopy(m,Ax,1,Axy,1);
			daxpy(m,-1.0,y,1,Axy,1);

			//funVal(iterStep)=Axy'*Axy/2;
			funVal[iterStep] = ddot(m,Axy,1,Axy,1)/2.0;
			for(i = 0 ; i<n ; i++)
				funVal[iterStep] += lambda * dnrm(k,&x[i],n,q);

			fVal = fValp = funVal[iterStep];
			if(iterStep != 0)
				fValp = funVal[iterStep-1];

			if(terminationCondition(&opt,fVal,fValp,x,dnrm2(n*k,xxp,1),dnrm2(n*k,xp,1),iterStep))
				break;
		}
	}
	ret.funVal = funVal;
	ret.errorCode = 0;
	ret.totIter = iterStep;
	strcpy(ret.type,"mtLeastRm1");
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;
}




FUNVAL mtLeastRm2(double* A, double* y,double z,OPTS opt){
	int m,n,k,i, iterStep;
	int *ind;
	double q, q_bar, lambda, lambda_max, alphap, alpha, L, beta ;
	double l_sum, r_sum, *funVal, *ValueL, fVal, fValp;
	double *tt, *ATy, *x, *Ax, *xp, *xxp, *t, *tp, *s, *As, *Axp, *ATAs, 
		*s_t, *u, *v, *Av, *g;
	double *Axy;
	FUNVAL ret;


	m = opt.m;
	n = opt.n;
	if(!hasGot(&opt,"ind"))
		errorOut("In mtLeastR, 'ind' should be specified.");
	else{
		ind = opt.ind;
		k = opt.ind_size - 1;
		if(ind[k] != (m))
			errorOut("Check opt.ind and its size.");
	}

	tt = (double *) malloc( sizeof(double) * n);
	ATy = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc( sizeof(double) * n*k);
	Ax = (double *) malloc( sizeof(double) * m);
	xp = (double *) malloc( sizeof(double) * n*k);
	Axp = (double *) malloc( sizeof(double) * m);
	xxp = (double *) malloc( sizeof(double) * n*k);
	t = (double *) malloc( sizeof(double) * n);
	tp = (double *) malloc( sizeof(double) * n);
	s = (double *) malloc( sizeof(double) * n*k);
	s_t = (double *) malloc( sizeof(double) * n);
	As = (double *) malloc( sizeof(double) * m);
	ATAs = (double *) malloc( sizeof(double) * n*k);
	g = (double *) malloc( sizeof(double) * n*k);
	u = (double *) malloc( sizeof(double) * n*k);
	v = (double *) malloc( sizeof(double) * n*k);
	Av = (double *) malloc( sizeof(double) * m);
	ValueL = (double *) malloc( sizeof(double) * opt.maxIter);
	Axy = (double *) malloc( sizeof(double) * m);
	funVal = (double *) malloc( sizeof(double) * opt.maxIter);

	if(!hasGot(&opt,"q"))
		opt.q = q = 2;
	else if( (q = opt.q) < 1)
		errorOut("'q' should be larger than 1.");
	
	initNormalization(&opt,A);

	/******* Starting point initialization *******/
/*
ATy=zeros(n, k);
% compute AT y
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    tt =A(ind_i,:)'*y(ind_i,1); / normalized
    ATy(:,i)= tt;
end
*/
	for(i=0; i < k; i++){
		normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&y[ind[i]],tt,ind[i],0,&opt);
		dcopy(n,tt,1,&ATy[i*n],1);
	}

	// process the regularization parameter
	if(opt.rFlag == 0)
		lambda = z;
	else{
		if( z<0.0 || z>1.0)
			errorOut("Z is not in range (0,1)");
		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATy[i],n,q_bar));
		lambda = z * lambda_max;
	}

	if( opt.init == 2)
		memset(x,0,sizeof(double)*n*k);
	else{
		if( hasGot(&opt,"x0"))
			x = opt.x0;
		else
			x = ATy;
	}

	//memset(Ax,0,sizeof(double)*m); //this is a waste of time actually..
	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0,&opt);

	if(opt.init == 0) {
		// TO-DO
	}

	/*........ Method specific code starts here .........*/

	/*************** Reformulated problem + Nemirovski's line search ******/
	if(opt.mFlag == 1 && opt.lFlag == 0 && opt.q == 2){
		L = 1.0;
		//xp=x; Axp=Ax; xxp=zeros(n,k);
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);
		memset(xxp,0,sizeof(double)*n*k);

		alphap=0; alpha=1;
		
		//t=sqrt(sum(x.^2, 2)); tp=t;
		for(i=0; i<n ; i++)
			t[i] = dnrm2(k,&x[i],n);
		dcopy(n,t,1,tp,1);

		for(iterStep = 0; iterStep < opt.maxIter ; iterStep ++ ){
			beta=(alphap-1)/alpha;    
			//s=x + beta* xxp; s_t= t + beta * (t-tp);
			dcopy(n*k,x,1,s,1);
			daxpy(n*k,beta,xxp,1,s,1);
			dcopy(n,t,1,s_t,1);
			dscal(n,1.0 + beta, s_t,1);
			daxpy(n,-beta,tp,1,s_t,1);

			//As=Ax + beta* (Ax-Axp);
			dcopy(m,Ax,1,As,1);
			dscal(m,1.0 + beta,As,1);
			daxpy(m, - beta,Axp,1,As,1);

			// COmputing ATAs
			for(i=0; i < k; i++){
				normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&As[ind[i]],tt,ind[i],0,&opt);
				dcopy(n,tt,1,&ATAs[i*n],1);
			}

			//g=ATAs-ATy;
			dcopy(n*k,ATAs,1,g,1);
			daxpy(n*k,-1.0,ATy,1,g,1);

			//xp=x;    Axp=Ax; tp=t;
			dcopy(n*k,x,1,xp,1);
			dcopy(m,Ax,1,Axp,1);
			dcopy(n,t,1,tp,1);

			for(;;){

				//u=s-g/L;
				dcopy(n*k,s,1,u,1);
				daxpy(n*k,-1.0/L,g,1,u,1);

				//v= s_t - lambda / L;
				dcopy(n,s_t,1,v,1);
				daxpy(n,-lambda/L,oneVector,1,v,1);

				ep21R(x,t,u,v,n,k);

				// v= x - s;
				dcopy(n*k,x,1,v,1);
				daxpy(n*k,-1.0,s,1,v,1);

				// computing Ax
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0,&opt);

				//Av=Ax -As;
				dcopy(m,Ax,1,Av,1);
				daxpy(m,-1.0,As,1,Av,1);

				//r_sum=norm(v,'fro')^2 + norm(t-s_t)^2; l_sum=norm(Av,'fro')^2;
				r_sum = ddot(n*k,v,1,v,1) -2.0*ddot(n,t,1,s_t,1) + ddot(n,t,1,t,1) + ddot(n,s_t,1,s_t,1);
				l_sum = ddot(m,Av,1,Av,1);

				if(l_sum <= r_sum * L)
					break;
				else
					L=maxof2(2.0*L, l_sum/r_sum);
			}
			ValueL[iterStep]=L;

			alphap = alpha; 
			alpha  = (1.0l+ sqrt(4.0l*alpha*alpha +1.0l))/2.0l;

			//xxp=x-xp;   Axy=Ax-y;
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1.0,xp,1,xxp,1);
			dcopy(m,Ax,1,Axy,1);
			daxpy(m,-1.0,y,1,Axy,1);
			//funVal(iterStep)=Axy'*Axy/2 +  lambda * sum(t);
			fVal = fValp = funVal[iterStep] = lambda*ddot(n,t,1,oneVector,1) + ddot(m,Axy,1,Axy,1)/2.0;

			if(iterStep != 0)
				fValp = funVal[iterStep - 1];

			if(terminationCondition(&opt,fVal,fValp,x,dnrm2(n*k,xxp,1),dnrm2(n*k,xp,1),iterStep))
				break;
		}
	}
	ret.errorCode=0;
	ret.funVal = funVal;
	ret.totIter = iterStep;
	strcpy(ret.type,"mtLeastRm2");
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;
}



FUNVAL mtLeastRm3(double* A, double* y,double z,OPTS opt){

	FUNVAL ret;
	int m,n,k,i, iterStep;
	int *ind;
	double q,q_bar, lambda, lambda_max, gamma, alphap, alpha, beta, L, tao;
	double fVal, fValp, *funVal, *ValueL, l_sum, r_sum;
	double *tt, *ATy, *x, *Ax, *s, *s_t, *As, *tp, *t, *xp, *xxp, *Axp, *ATAs, *ATAx, 
		*ATAxp, *Axy, *tnew, *xnew, *Axnew, *v, *v_t, *Av, *u, *g ;

	m = opt.m;
	n = opt.n;
	if(!hasGot(&opt,"ind"))
		errorOut("In mtLeastR, 'ind' should be specified.");
	else{
		ind = opt.ind;
		k = opt.ind_size - 1;
		if(ind[k] != (m))
			errorOut("Check opt.ind and its size.");
	}

	tt = (double *) malloc( sizeof(double) * n);
	ATy = (double *) malloc( sizeof(double) * n*k);
	x = (double *) malloc( sizeof(double) * n*k);
	Ax = (double *) malloc( sizeof(double) * m);
	xp = (double *) malloc( sizeof(double) * n*k);
	Axp = (double *) malloc( sizeof(double) * m);
	xxp = (double *) malloc( sizeof(double) * n*k);
	t = (double *) malloc( sizeof(double) * n);
	tp = (double *) malloc( sizeof(double) * n);
	As = (double *) malloc( sizeof(double) * m);
	ATAs = (double *) malloc( sizeof(double) * n*k);
	ATAx = (double *) malloc( sizeof(double) * n*k);
	ATAxp = (double *) malloc( sizeof(double) * n*k);
	s = (double *) malloc( sizeof(double) * n*k);
	s_t = (double *) malloc( sizeof(double) * n);
	g = (double *) malloc( sizeof(double) * n*k);
	u = (double *) malloc( sizeof(double) * n*k);
	v = (double *) malloc( sizeof(double) * n*k);
	xnew = (double *) malloc( sizeof(double) * n*k);
	v_t = (double *) malloc( sizeof(double) * n);
	tnew = (double *) malloc( sizeof(double) * n);
	Axnew = (double *) malloc( sizeof(double) * m);
	Av = (double *) malloc( sizeof(double) * m);
	ValueL = (double *) malloc( sizeof(double) * opt.maxIter);
	Axy = (double *) malloc( sizeof(double) * m);
	funVal = (double *) malloc( sizeof(double) * opt.maxIter);


	if(!hasGot(&opt,"q"))
		opt.q = q = 2;
	else if( (q = opt.q) < 1)
		errorOut("'q' should be larger than 1.");
	
	initNormalization(&opt,A);

	/******* Starting point initialization *******/
/*
ATy=zeros(n, k);
% compute AT y
for i=1:k
    ind_i=(ind(i)+1):ind(i+1);     % indices for the i-th group
    tt =A(ind_i,:)'*y(ind_i,1); / normalized
    ATy(:,i)= tt;
end
*/
	for(i=0; i < k; i++){
		normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&y[ind[i]],tt,ind[i],0,&opt);
		dcopy(n,tt,1,&ATy[i*n],1);
	}

	// process the regularization parameter
	if(opt.rFlag == 0)
		lambda = z;
	else{
		if( z<0.0 || z>1.0)
			errorOut("Z is not in range (0,1)");
		if(q==1)
			q_bar = DBL_MAX;
		else if ( q> 1.0e6 )
			q_bar = 1.0;
		else
			q_bar = q/(q-1);
		lambda_max = 0.0;
		for(i=0; i<n ; i++)
			lambda_max = maxof2(lambda_max ,dnrm(k,&ATy[i],n,q_bar));
		lambda = z * lambda_max;
	}

	if( opt.init == 2)
		memset(x,0,sizeof(double)*n*k);
	else{
		if( hasGot(&opt,"x0"))
			x = opt.x0;
		else
			x = ATy;
	}

	memset(Ax,0,sizeof(double)*m);
	for(i=0; i < k; i++)
		normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&x[i*n],&Ax[ind[i]],ind[i],0,&opt);

	if(opt.init == 0) {
		// TO-DO
	}

	/*........ Method specific code starts here .........*/

	/********* Reformulated problem + adaptive Line Search *****/
	if (opt.mFlag==1 && opt.lFlag==1 && opt.q==2){
		L = 1.0;
		gamma = 1.0;
		//xp=x; Axp=Ax;
		dcopy(n*k,x,1,xp,1);
		dcopy(m,Ax,1,Axp,1);

		memset(xxp,0,sizeof(double)*n*k);

		//t=sqrt(sum(x.^2, 2)); tp=t;
		for(i=0; i<n ; i++)
			t[i] = dnrm2(k,&x[i],n);
		dcopy(n,t,1,tp,1);

		//As= Ax; 
		dcopy(m,Ax,1,As,1);

		// COmputing ATAs
		for(i=0; i < k; i++){
			normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&As[ind[i]],tt,ind[i],0,&opt);
			dcopy(n,tt,1,&ATAs[i*n],1);
		}

		//ATAx=ATAs;
		dcopy(n*k,ATAs,1,ATAx,1);

		for(iterStep = 0; iterStep < opt.maxIter ; iterStep ++ ){
			//ATAxp=ATAx;
			dcopy(n*k,ATAx,1,ATAxp,1);

			// As = Ax;
			dcopy(m,Ax,1,As,1);
			// Computing ATAs
			for(i=0; i < k; i++){
				normalizedmv_lowlevel('T',ind[i+1]-ind[i],n,&A[ind[i]],m,&As[ind[i]],tt,ind[i],0,&opt);
				dcopy(n,tt,1,&ATAs[i*n],1);
			}

			//ATAx=ATAs;
			dcopy(n*k,ATAs,1,ATAx,1);

			// Line Search for L begins
			for(;;){
				if(iterStep != 0){
	                alpha= (-gamma+ sqrt(gamma*gamma + 4.0* L * gamma)) / (2.0*L);
		            beta= (gamma - gamma* alphap) / (alphap * gamma + alphap* L * alpha);
					//        % beta is the coefficient for generating search point s

					//s=x + beta* xxp;   s_t= t + beta * (t -tp);
					dcopy(n*k,x,1,s,1);
					daxpy(n*k,beta,xxp,1,s,1);
					dcopy(n,t,1,s_t,1);
					dscal(n,1.0 + beta,s_t,1);
					daxpy(n,-beta,tp,1,s_t,1);

					//As=Ax + beta* (Ax-Axp);
					dcopy(m,Ax,1,As,1);
					dscal(m,1.0+beta,As,1);
					daxpy(m,-beta,Axp,1,As,1);

					//ATAs=ATAx + beta * (ATAx- ATAxp);
					dcopy(n*k,ATAx,1,ATAs,1);
					dscal(n*k,1.0 + beta,ATAs,1);
					daxpy(n*k,-beta,ATAxp,1,ATAs,1);
				}
				else{
	                alpha= (-1.0l+ sqrt(5.0l)) / 2.0l;
		            beta=0.0l; 
			        //s=x; s_t=t; As=Ax; ATAs=ATAx;
					dcopy(n*k,x,1,s,1);
					dcopy(n,t,1,s_t,1);
					dcopy(m,Ax,1,As,1);
					dcopy(n*k,ATAx,1,ATAs,1);
				}

				//g=ATAs-ATy;
				dcopy(n*k,ATAs,1,g,1);
				daxpy(n*k,-1.0,ATy,1,g,1);

				//u=s-g / L;
				dcopy(n*k,s,1,u,1);
				daxpy(n*k,-1.0/L,g,1,u,1);

				//v= s_t - lambda / L;
				dcopy(n,s_t,1,v,1);
				daxpy(n,-lambda/L,oneVector,1,v,1);

				ep21R(xnew,tnew,u,v,n,k);

				//v=xnew-s;
				dcopy(n*k,xnew,1,v,1);
				daxpy(n*k,-1.0,s,1,v,1);

				//v_t=tnew-s_t;
				dcopy(n,tnew,1,v_t,1);
				daxpy(n,-1.0,s_t,1,v_t,1);

				//memset(Axnew,0,sizeof(double)*m); // waste of time
				for(i=0; i < k; i++)
					normalizedmv_lowlevel('N',ind[i+1]-ind[i],n,&A[ind[i]],m,&xnew[i*n],&Axnew[ind[i]],ind[i],0,&opt);
				

				//Av=Axnew -As;
				dcopy(m,Axnew,1,Av,1);
				daxpy(m,-1.0,As,1,Av,1);
				//r_sum=norm(v,'fro')^2 + v_t'*v_t; l_sum=norm(Av,'fro')^2;
				r_sum = ddot(n*k,v,1,v,1) + ddot(n,v_t,1,v_t,1);
				l_sum = ddot(m,Av,1,Av,1);	
				if(l_sum <= r_sum * L)
					break;
				else
					L=max(2*L, l_sum/r_sum);
			}
			gamma=L* alpha* alpha;    
			alphap=alpha;

			ValueL[iterStep]=L;

			tao=L * r_sum / l_sum;
			if (tao >= 5.0)
				L=L*0.8;

	        //xp=x;    x=xnew; xxp=x-xp;  
		    //Axp=Ax;  Ax=Axnew;

			dcopy(n*k,x,1,xp,1);
			dcopy(n*k,xnew,1,x,1);
			dcopy(n*k,x,1,xxp,1);
			daxpy(n*k,-1.0,xp,1,xxp,1);
			dcopy(m,Ax,1,Axp,1);
			dcopy(m,Axnew,1,Ax,1);

			//tp=t; t=tnew;
			dcopy(n,t,1,tp,1);
			dcopy(n,tnew,1,t,1);

			//Axy=Ax-y;
			dcopy(m,Ax,1,Axy,1);
			daxpy(m,-1.0,y,1,Axy,1);

			//funVal(iterStep)=Axy'*Axy/2 + lambda * sum(t);
			fValp = fVal = funVal[iterStep] = lambda*ddot(n,t,1,oneVector,1) + ddot(m,Axy,1,Axy,1)/2.0;
			if(iterStep != 0 )
				fValp = funVal[iterStep -1];

			if(terminationCondition(&opt,fVal,fValp,x,dnrm2(n*k,xxp,1),dnrm2(n*k,xp,1),iterStep))
				break;
		}
	}
	ret.errorCode=0;
	ret.funVal = funVal;
	ret.totIter = iterStep;
	strcpy(ret.type,"mtLeastRm2");
	ret.ValueL = ValueL;
	ret.x = x;
	return ret;
}